4_A549_27plex_analysis

Down stream analysis of the dataset ‘4_A549_27plex’

setup

knitr::opts_chunk$set(dev = "png",dpi = 80)
library(ggforce)
library(tidyverse)

dcolor <- ggsci::scale_color_d3('category20')
dfill <- ggsci::scale_fill_d3('category20')
ccolor <- viridis::scale_color_viridis()
style <- function(x,color,method='UMAP',axis=c(1,2),Theme=NULL,
                  axis.text=FALSE,axis.title=FALSE,coord_fix =TRUE,
                  palette=TRUE,title=NULL,
                  legend=TRUE,scale_color_log10=FALSE,
                  guide_pointSize=2.5,...) {
  axis <- paste0(method,axis)
  if(is.null(Theme)) Theme <- theme_bw
  g <- ggplot() + Theme()
  g <- g + geom_point(aes_(x=as.name(axis[1]),y=as.name(axis[2]),color=as.name(color)),x,...)
  
  if(coord_fix) {
    g <- g + coord_fixed()
  }
  if(!axis.text) {
    g <- g + theme(axis.text = element_blank(),axis.ticks = element_blank())
  }
  if(!axis.title) {
    g <- g + theme(axis.title = element_blank())
  }
  if(is.null(title)) {
    g <- g + ggtitle(method) 
  }else if(!is.na(title)) {
    g <- g + ggtitle(title) 
  }
  if(!legend) g <- g + theme(legend.position = 'none')
  
  if(scale_type(as.matrix(x[,color]))=='discrete' & !is.null(guide_pointSize)) {
    g <- g + guides(color = guide_legend(override.aes = list(size=guide_pointSize)))
  }
  
  if(palette) {
    if(scale_type(as.matrix(x[,color]))=='discrete') {
      g <- g + ggsci::scale_color_d3('category20')
    }else{
      if(scale_color_log10){
        g <- g + viridis::scale_color_viridis(trans='log10')
      }else{
        g <- g + viridis::scale_color_viridis()
      }
    }
  }
  return(g)
}
tib2df <- function(tib,RowNames=1) {
  rn <- pull(tib,all_of(RowNames))
  df <- tib %>%
    select(-RowNames) %>%
    as.data.frame()
  rownames(df) <- rn
  return(df)
}

load data

data <- tibble(
  path = list.files("4_A549_27plex/Raw/",full.names = T),
  sampleName = c('0h-r1', '0h-r2', '0p5h-r1', '0p5h-r2', '24h-r1', '24h-r2', '48h-r1', '48h-r2')
  ) %>%
  separate(sampleName,into=c('treat','rep'),remove = F) %>%
  mutate(
    metric = map(path,~{
      .x %>%
        read_csv(show_col_types=F) %>% 
        unite(cell,tileID,cellID,sep = '-',remove = F)
      }),
    IF = map(metric,~{.x[,-(2:8)]}),
    meta = map2(metric,IF,~{
      .x[,1:8] %>%
        unite(cell,tileID,cellID,sep = '-',remove = F) %>%
        mutate(totalExp = rowSums(.y[,-1]))
      }),
    ) %>%
  select(-path,-metric)

Quality check & control

data %>%
  select(-IF) %>% unnest(meta) %>%
  ggplot(aes(x,y,color=treat)) +
    theme_void() +
    geom_point(size=.1) + dcolor +
    facet_grid(rep~treat) + coord_fixed()

data %>%
  select(sampleName,meta) %>% unnest(meta) %>%
  ggplot(aes(sampleName,totalExp)) + geom_violin() + geom_boxplot()

dataf <- data %>%
  mutate(
    idx = map(meta,~{.x$totalExp > 1.3}),
    IF = map2(IF,idx,~{.x[.y,]}),
    meta = map2(meta,idx,~{.x[.y,]})
    ) %>%
  mutate(
    idx2 = map(IF,~{rowSums(.x[,-1]==0) == 0}),
    IF = map2(IF,idx2,~{.x[.y,]}),
    meta = map2(meta,idx2,~{.x[.y,]}),
    nCell = map(idx2,sum) %>% unlist()
    ) %>%
  select(-idx,-idx2)

dataf %>%
  select(sampleName,meta) %>% unnest(meta) %>%
  ggplot(aes(sampleName,totalExp)) + 
  geom_violin() + geom_boxplot()

dataf %>%  
  ggplot(aes(sampleName,nCell)) + 
  theme_classic() +
  geom_bar(stat='identity') +
  scale_x_discrete(limits=rev) +
  theme(axis.title.y = element_blank()) +  
  coord_flip() + theme(aspect.ratio = .35) 

tmp <- dataf %>% select(sampleName,IF) %>% unnest(IF)
idx <- with(tmp,model.matrix(~0+sampleName))
colnames(idx) <- sub("sampleName","",colnames(idx))
avgmat <- tmp[,-(1:2)] %>% as.matrix() %>% t %>% {t(.%*%idx)/colSums(idx)} %>% scale %>% asinh %>% t 
pheatmap::pheatmap(avgmat,cluster_cols = F,cellwidth = 10,color=viridis::viridis(10),show_rownames = T)

Meta <- dataf %>%
  select(sampleName,treat,rep,meta) %>% unnest(meta) %>%
  unite(Cell,treat,rep,cell,sep = '-',remove = F)
X <- bind_rows(dataf$IF) %>% select(-cell) %>% 
   scale() %>% asinh()
rownames(X) <- Meta$Cell

p <- prcomp(X)
screeplot(p,type='l',n=30)

total <- bind_rows(dataf$meta) %>% pull(totalExp) 
cor(p$x,total) %>%
  plot(xlab = 'PC',ylab = 'CorCoef_vsTotalExp')

um <- uwot::umap(p$x[,2:6],
                 metric = 'cosine',
                 n_neighbors = 30,
                 min_dist = .1) %>%
  as_tibble(rownames = 'Cell') %>%
  dplyr::rename(UMAP1=2,UMAP2=3)

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
Using compatibility `.name_repair`.

Metaum <- inner_join(Meta,um,by = 'Cell')

style(Metaum,'treat',title = NULL,size = .1,alpha = .6)

Warning: `aes_()` was deprecated in ggplot2 3.0.0.
Please use tidy evaluation idioms with `aes()`

style(Metaum,'sampleName',title = NA,size = .1,alpha = .6)

style(Metaum,'totalExp',title = NA,size = .1,alpha = .6)

Pseudo time

Using the datasets of 0h and 0.5h treatment

Meta2 <- Meta %>% filter(treat %in% c('0h','0p5h'))
X2 <- dataf %>%
  filter(treat %in% c('0h','0p5h')) %>%
  pull(IF) %>% bind_rows() %>% 
  select(-cell) %>% 
  scale() %>% asinh()
rownames(X2) <- Meta2$Cell

p2 <- prcomp(X2)
screeplot(p2,type='l',n=30)

total2 <- Meta2$totalExp
cor(p2$x,total2) %>%
  plot(xlab = 'PC',ylab = 'CorCoef_vsTotalExp')

um2 <- uwot::umap(p2$x[,2:6],
                 metric = 'cosine',
                 n_neighbors = 30,
                 min_dist = .1) %>%
  as_tibble(rownames = 'Cell') %>%
  dplyr::rename(UMAP1=2,UMAP2=3)
Metaum2 <- inner_join(Meta2,um2,by = 'Cell')
style(Metaum2,'treat',title = NULL,size = .1,alpha = .6)

ph2 <- phateR::phate(p2$x[,2:6],
                     mds.dist.method = 'cosine',
                     knn = 20)

Calculating PHATE...
  Running PHATE on 19597 observations and 5 variables.
  Calculating graph and diffusion operator...
    Calculating KNN search...
    Calculated KNN search in 2.02 seconds.
    Calculating affinities...
    Calculated affinities in 0.10 seconds.
  Calculated graph and diffusion operator in 2.14 seconds.
  Calculating landmark operator...
    Calculating SVD...
    Calculated SVD in 1.36 seconds.
    Calculating KMeans...
    Calculated KMeans in 2.33 seconds.
  Calculated landmark operator in 4.20 seconds.
  Calculating optimal t...
    Automatically selected t = 26
  Calculated optimal t in 2.87 seconds.
  Calculating diffusion potential...
  Calculated diffusion potential in 0.58 seconds.
  Calculating metric MDS...
  Calculated metric MDS in 5.23 seconds.
Calculated PHATE in 15.02 seconds.

Metaum2 <- Metaum2 %>%
  mutate(PHATE1 = ph2$embedding[,1],PHATE2 = ph2$embedding[,2])
style(Metaum2,'treat','PHATE',size=.1)

style(Metaum2,'rep','PHATE',size=.1,palette = F)

set.seed(1); km <- kmeans(ph2$embedding,centers=5)
Metaum2 %>% 
  mutate(km=as.factor(km$cluster)) %>%
  style('km','PHATE',size=.1)

library(slingshot)

Loading required package: princurve
Loading required package: TrajectoryUtils
Loading required package: SingleCellExperiment
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following object is masked _by_ ‘.GlobalEnv’:

    count

The following object is masked from ‘package:dplyr’:

    count


Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps,
    colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates,
    colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse,
    rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans,
    rowWeightedMedians, rowWeightedSds, rowWeightedVars

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians


Attaching package: ‘SummarizedExperiment’

The following object is masked from ‘package:Seurat’:

    Assays

The following object is masked from ‘package:SeuratObject’:

    Assays

sce <- getLineages(ph2$embedding,as.character(km$cluster),start.clus = '5') %>%
  getCurves() %>% slingPseudotime() %>%
  as_tibble(rownames='Cell') %>% mutate(Lineage1 = Lineage1/max(Lineage1)) %>%
  rename(PseudoTime = Lineage1)

Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.

Metaum2 <- left_join(Metaum2,sce,by='Cell')
Metaum2 %>% style("PseudoTime","PHATE",size=.2) + theme_void()

Metaum2 %>% style("PseudoTime","UMAP",size=.2) + theme_void()

Metaum2 %>%
  ggplot(aes(PseudoTime,fill=treat)) + geom_histogram(position = 'dodge') + theme_classic() + dfill + coord_fixed(ratio=.0003) +
  theme(axis.ticks = element_blank(),panel.grid = element_blank())

LOESS <- dataf %>%
  filter(treat %in% c("0h","0p5h")) %>%
  select(sampleName,IF) %>% unnest(IF) %>%
  unite(Cell,sampleName,cell,sep='-') %>%
  left_join(select(Metaum2,Cell,PseudoTime),.) %>%
  gather(Prot,ExpLevel,-(1:2)) %>%
  group_by(Prot) %>% nest() %>%
  mutate(LOESS = map(data,~loess(ExpLevel~PseudoTime,.x,span = 0.2)),
         PRED = map(LOESS,~predict(.x,newdata=seq(0,1,by=0.01))),
         newdata = map(PRED,~tibble(PseudoTime=seq(0,1,by=0.01),PRED=.x)) )

Joining with `by = join_by(Cell)`

loessPl <- function(x,ratio = .35){
  ggplot(x,aes(PseudoTime,scaled,color=Prot)) + 
  geom_line() + theme_bw() +
  theme(legend.position = 'none') +
  theme(axis.ticks.y = element_blank(),
        strip.background = element_rect(fill=NA,color=NA), 
        panel.grid.major.x = element_blank(),panel.grid.minor = element_blank()) + 
  coord_fixed(ratio=ratio) + labs(y = 'scaled expression level')
}

LOESS_scaled <- LOESS %>% select(Prot,newdata) %>% 
  mutate(newdata = map(newdata,~mutate(.x,scaled = (PRED-min(PRED))/(max(PRED)-min(PRED)) ))) %>%
  unnest(newdata)  %>% ungroup()

usep <- c("p4EBP1/2/3","pATM","pELF2A","pERK","pFAK","pp38","pSMAD1/5","pSMAD3","pSTAT3")
LOESS_scaled %>% filter(Prot %in% usep) %>% loessPl()

LOESS_scaled %>% filter(Prot %in% usep) %>% loessPl() + facet_wrap(~Prot)

LOESS_scaled %>% filter(!Prot %in% usep) %>% loessPl() + facet_wrap(~Prot)

LOESS_scaled %>% filter(!Prot %in% usep) %>% loessPl()

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] slingshot_2.6.0             TrajectoryUtils_1.6.0       SingleCellExperiment_1.20.1 SummarizedExperiment_1.28.0 Biobase_2.58.0              MatrixGenerics_1.10.0       matrixStats_1.0.0          
 [8] princurve_2.1.6             lubridate_1.9.2             forcats_1.0.0               stringr_1.5.0               dplyr_1.1.2                 purrr_1.0.1                 readr_2.1.4                
[15] tidyr_1.3.0                 tibble_3.2.1                tidyverse_2.0.0             ggforce_0.4.1               ggplot2_3.4.2               Seurat_4.3.0.1              SeuratObject_4.1.3         
[22] sp_2.0-0                    GenomicRanges_1.50.2        GenomeInfoDb_1.34.9         IRanges_2.32.0              S4Vectors_0.36.2            BiocGenerics_0.44.0        

loaded via a namespace (and not attached):
  [1] utf8_1.2.3                spatstat.explore_3.2-1    reticulate_1.30           tidyselect_1.2.0          RSQLite_2.3.1             AnnotationDbi_1.60.2      htmlwidgets_1.6.2        
  [8] grid_4.2.2                Rtsne_0.16                munsell_0.5.0             codetools_0.2-19          ica_1.0-3                 future_1.33.0             miniUI_0.1.1.1           
 [15] withr_2.5.0               spatstat.random_3.1-5     colorspace_2.1-0          progressr_0.13.0          knitr_1.43                rstudioapi_0.15.0         ROCR_1.0-11              
 [22] tensor_1.5                listenv_0.9.0             labeling_0.4.2            GenomeInfoDbData_1.2.9    polyclip_1.10-4           bit64_4.0.5               farver_2.1.1             
 [29] pheatmap_1.0.12           rprojroot_2.0.3           parallelly_1.36.0         vctrs_0.6.3               generics_0.1.3            xfun_0.39                 timechange_0.2.0         
 [36] R6_2.5.1                  phateR_1.0.7              bitops_1.0-7              spatstat.utils_3.0-3      cachem_1.0.8              DelayedArray_0.24.0       promises_1.2.0.1         
 [43] scales_1.2.1              vroom_1.6.3               gtable_0.3.3              globals_0.16.2            goftest_1.2-3             rlang_1.1.1               pamr_1.56.1              
 [50] splines_4.2.2             lazyeval_0.2.2            spatstat.geom_3.2-4       yaml_2.3.7                reshape2_1.4.4            abind_1.4-5               stylo_0.7.4              
 [57] httpuv_1.6.11             tools_4.2.2               tcltk_4.2.2               ellipsis_0.3.2            jquerylib_0.1.4           RColorBrewer_1.1-3        ggdendro_0.1.23          
 [64] proxy_0.4-27              ggridges_0.5.4            Rcpp_1.0.11               plyr_1.8.8                sparseMatrixStats_1.10.0  zlibbioc_1.44.0           RCurl_1.98-1.12          
 [71] deldir_1.0-9              pbapply_1.7-2             viridis_0.6.4             cowplot_1.1.1             zoo_1.8-12                ggrepel_0.9.3             cluster_2.1.4            
 [78] here_1.0.1                magrittr_2.0.3            data.table_1.14.8         scattermore_1.2           lmtest_0.9-40             RANN_2.6.1                fitdistrplus_1.1-11      
 [85] hms_1.1.3                 patchwork_1.1.2           mime_0.12                 evaluate_0.21             xtable_1.8-4              gridExtra_2.3             compiler_4.2.2           
 [92] KernSmooth_2.23-22        crayon_1.5.2              htmltools_0.5.5           later_1.3.1               tzdb_0.4.0                DBI_1.1.3                 tweenr_2.0.2             
 [99] MASS_7.3-60               Matrix_1.5-4.1            cli_3.6.1                 parallel_4.2.2            igraph_1.5.0.1            pkgconfig_2.0.3           plotly_4.10.2            
[106] spatstat.sparse_3.0-2     bslib_0.5.0               XVector_0.38.0            digest_0.6.33             sctransform_0.3.5         RcppAnnoy_0.0.21          tsne_0.1-3.1             
[113] spatstat.data_3.0-1       Biostrings_2.66.0         rmarkdown_2.23            leiden_0.4.3              uwot_0.1.16               DelayedMatrixStats_1.20.0 shiny_1.7.4.1            
[120] lifecycle_1.0.3           nlme_3.1-162              jsonlite_1.8.7            viridisLite_0.4.2         fansi_1.0.4               pillar_1.9.0              ggsci_3.0.0              
[127] lattice_0.21-8            KEGGREST_1.38.0           fastmap_1.1.1             httr_1.4.6                survival_3.5-5            glue_1.6.2                png_0.1-8                
[134] bit_4.0.5                 tcltk2_1.2-11             class_7.3-22              stringi_1.7.12            sass_0.4.7                blob_1.2.4                memoise_2.0.1            
[141] irlba_2.3.5.1             e1071_1.7-13              future.apply_1.11.0       ape_5.7-1